* Independently generated unit root outcome and unit root + white noise predictors DGP

* Set working directory
cd"/Users/alikagalwala/Library/CloudStorage/Dropbox/Whitten_Kagalwala/TS/JOP/Acceptance/KW_Replication/Multiple_Predictors/URDV-MixIVs"

* Run the file with the simulation program
do "URDV-MixIV-DGP.do"

// Monte Carlo Analyses
* Create a matrix to store results
matrix mc_summary_stats = J(7, 25, .)

matrix colnames mc_summary_stats = phi_gecm_t1 type1_x_gecm type1_lx_gecm phi_adl_t1 type1_x_adl11 type1_lx_adl11 ///
phi1_adl22_t1 type1_x_adl22 type1_lx_adl22 phi1_adl33_t1 type1_x_adl33 type1_lx_adl33 ///
type1_x1_gecm type1_x2_gecm type1_x1_adl11 type1_x2_adl11 type1_x1_adl22 type1_x2_adl22 type1_x1_adl33 type1_x2_adl33 phi_gecm_power phi_adl_power phi1_adl22_power phi1_adl33_power time
 
   local j = 1
 foreach t in 25 50 75 100 150 500 1000 {
* Create a counter to store results in the matrix
parallel sim, expr(phi_gecm_t1 = r(phi_gecm_t1) type1_x_gecm = r(type1_x_gecm) type1_lx_gecm = r(type1_lx_gecm) ///
	phi_adl_t1 = r(phi_adl_t1) type1_x_adl11 = r(type1_x_adl11) type1_lx_adl11 = r(type1_lx_adl11) ///
	phi1_adl22_t1 = r(phi1_adl22_t1) type1_x_adl22 = r(type1_x_adl22) type1_lx_adl22 = r(type1_lx_adl22) ///
	phi1_adl33_t1 = r(phi1_adl33_t1) type1_x_adl33 = r(type1_x_adl33)  type1_lx_adl33 = r(type1_lx_adl33) ///
	type1_x1_gecm = r(type1_x1_gecm) type1_x2_gecm = r(type1_x2_gecm) type1_x1_adl11 = r(type1_x1_adl11) ///
	type1_x2_adl11 = r(type1_x2_adl11) type1_x1_adl22 = r(type1_x1_adl22)  type1_x2_adl22 = r(type1_x2_adl22) type1_x1_adl33 = r(type1_x1_adl33) type1_x2_adl33 = r(type1_x2_adl33) phi_gecm_power = r(phi_gecm_power) phi_adl_power = r(phi_adl_power) ///
					phi1_adl22_power = r(phi1_adl22_power) phi1_adl33_power = r(phi1_adl33_power)) ///
	reps(1000) noisily seeds(12344 12344 12344 12344 12344 12344 12344 12344): NG, obs(`t')
			
					qui su phi_gecm_t1
					matrix mc_summary_stats[`j',1] = r(mean)
					qui su type1_x_gecm
					matrix mc_summary_stats[`j',2] = r(mean)
					qui su type1_lx_gecm
					matrix mc_summary_stats[`j',3] = r(mean)
					qui su phi_adl_t1
					matrix mc_summary_stats[`j',4] = r(mean)
					qui su type1_x_adl11
					matrix mc_summary_stats[`j',5] = r(mean)
					qui su type1_lx_adl11
					matrix mc_summary_stats[`j',6] = r(mean)
					qui su phi1_adl22_t1
					matrix mc_summary_stats[`j',7] = r(mean)
					qui su type1_x_adl22
					matrix mc_summary_stats[`j',8] = r(mean)
					qui su type1_lx_adl22
					matrix mc_summary_stats[`j',9] = r(mean)
					qui su phi1_adl33_t1
					matrix mc_summary_stats[`j',10] = r(mean)
					qui su type1_x_adl33
					matrix mc_summary_stats[`j',11] = r(mean)
					qui su type1_lx_adl33
					matrix mc_summary_stats[`j',12] = r(mean)
					qui su type1_x1_gecm
					matrix mc_summary_stats[`j',13] = r(mean)
					qui su type1_x2_gecm
					matrix mc_summary_stats[`j',14] = r(mean)
					qui su type1_x1_adl11
					matrix mc_summary_stats[`j',15] = r(mean)
					qui su type1_x2_adl11
					matrix mc_summary_stats[`j',16] = r(mean)
					qui su type1_x1_adl22
					matrix mc_summary_stats[`j',17] = r(mean)
					qui su type1_x2_adl22
					matrix mc_summary_stats[`j',18] = r(mean)
					qui su type1_x1_adl33
					matrix mc_summary_stats[`j',19] = r(mean)
					qui su type1_x2_adl33
					matrix mc_summary_stats[`j',20] = r(mean)

					* Power
					su phi_gecm_power
					matrix mc_summary_stats[`j', 21] = r(mean)
					su phi_adl_power
					matrix mc_summary_stats[`j', 22] = r(mean)
					su phi1_adl22_power
					matrix mc_summary_stats[`j', 23] = r(mean)
					su phi1_adl33_power
					matrix mc_summary_stats[`j', 24] = r(mean)
					
					matrix mc_summary_stats[`j', 25] = `t'
					local j = `j' + 1
 }
 
// Store MC matrix results as a dataframe
clear
set obs 7
svmat double mc_summary_stats, names(col)
save multiple_predictors_URDVMIXIVs.dta, replace

